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ABSTRACT 

We develop two different solar dynamo models to verify the hypothesis that a deep 
meridional flow can restrict the apperance of sunspots below 45 degrees, proposed by 
iNandv fc Choudhuril ([2002). In the first one, a single polytropic approximation for the 
density profile was taken, for both radiative and convective zo nes. In the second one, 
two p olytropes were used to distinguish between both zone s llPinzon fc Calvo-Mozo] 
120011) . The magnetic buoyancy mechanism proposed by iDikpati fc Charbonneaul 
(1999J) was chosen in both models. We, actually, have obtained that a deep merid- 
ional flow pus hes the maxima of toroidal m agnetic field toward the solar equator, but 
in contrast to INandv fc Choudhuril (2002), a second zone of maximal fields remains 
at the poles. The second model, although clos e ly resembling the solar standard mode l 
of Bahca lT Pinsonneault fc Wasserbual)1995|) : iBahcall. Pinsonneault fc Basul l)2001j) . 
gives solar cyles three times longer than observed. 

Key words: MHD - magnetic field - Sun: interior - Sun: magnetic fields. 



1 INTRODUCTION 

Despite its irregular appearance, the solar magnetic cycle 
evolves in a spatially and temporally well organized manner. 
Some sunspots associated phenomena, such as the 11 years 
cycle of alternating polarities (Hale law), allows to conclude 
that the solar magnetic cycle can be explained by a dynamo 
process involving the transformation of a poloidal magnetic 
field into a toroidal one, and the later regeneration of the 
poloidal field, but of opposed polarity to the initial, and 
so on. In the Babcock-Leighton approach iBabcocklll96lt 
iLeiehtonl 1969 s ! the dynamo operates in the following way: In 
a first stage, an initial dipolar field, with lines on meridional 
planes, is dragged in the east- west direction by the action of 
the solar differential rotation to form a toroidal field. This 
happens close to the base of the solar convection zone (SCZ), 
in a thin layer called tachocline, where helioseismology has 
discovered a substancial radial shear in the rotation pattern. 
In a second stage, tubes of toroidal flux emerge to the sur- 
face, since these tubes of intense magnetic field are less dense 
than their sourrondings (magnetic bouyancy force). While 
the force lines rises they are twisted by the Coriolis force to 
form the bipolar magnetic regions (BMR) associated with 
sunspots in the solar surface. In the last stage, the leader 
portion of each BMR (i.e, that is farther ahead in the direc- 
tion of solar rotation) migrates equatorward, while the fol- 
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lower portion migrates poleward, thanks to the meridional 
circulation. The general effect is that the force lines between 
the two BMRs lay on the poloidal direction, but oppossite 
to the original field. On the long run, they cancel out the 
original dipol ar field and give rise to a new one in the oppo- 
site direction Jch oudhuri. Schiisler fc Dik patill995l : IPurnevI 
Il995lll996l.ll997l: IDikpati fc Charbonneaulll999l) 

In the kinematic regime, where a fixed velocity field 
is assumed and one investigates just the evolution of the 
magnetic field due to that velocitiy field (without the 
back-reaction of the magnetic field onto the solar plasma) 
iChoudhuri200ol) , three main ingredients are used to resem- 
ble the mass transport in the Babcock-Leighton approach: 
differential rotation, meridional circulation and magnetic 
buoyancy. The first one can be suitably determined from 
helioseismology measurements, and resembles closely the 
radial shear at the tachocline and the latitudinal angu- 
lar velocity distribution in the convection zone. The sec- 
ond one is the main flux transport agent. It is responsible 
to lead the poloidal flux generated in the surface to the 
deeper layers where toroidal flux is regenerated. Obeser- 
vational evidence tell us that there is a poleward merid- 
ional flow at the surface with an average speed between 
10—20 m s" 1 , but the internal return flow remains unknown 
JCiles et al .119971) . The last ingredient (magnetic bouyancy) 
summarizes the result s of ris i ng magnetic tube simu l ations 
JP'Silva fc Choudhuril Il993t iFan. Fisher fc DeLucal Il993t 
Calieari, Moreno-Insertis fc Schiislerl Il995l Il998ft . where a 
buoyant magnetic force appears, thanks to the density gra- 
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dient between the inside and the outside of the tube. We will 
incorporate a simplified form for magnetic buoyancy which 
resembles an average number of such events. 

Most of the recent kinematic models are able to sucess- 
fully explain the phenomena like the reversal of the solar 
dipolar field in an 11-years cycle or the reversal of po- 
laritie s of sunspot pairs from one cycle to the next (Hale 
law) jDikpati fc Choudhurilll99i iDikpati fc Charbonneaul 
1999). Neverthless, they fail to predict the absence of 
sunspots at latitudes above 30 degrees (Sporer Law). 
Based on recent numerical simulations which suggest that 
a subadiabatic stratification, like the one in the so- 
lar radia tive zone, is able to resist the magnetic buoy- 
ancy jRempel. Schiiser fc T6tnll2000h . lNandv fc Choudhuril 
(2002) proposed a dynamo model which restrics for the first 
time the sunspots' appearance below 45 degrees. In this 
model, the meridional flow penetrates below the solar con- 
vection zone, where the toroidal magnetic field can't erupt 
to the surface, and is dragged in equatorward direction by 
the meridional circulation. All these models distinguish for 
the magnetic diffusivity and the differential rotation profiles 
between the radiative and convective zones, but use a merid- 
ional flow based on a single polytropic approximation to the 
density profile for both zones. 

This work i nvesti gates the hypothesis of 
iNandv fc Choudhuril i2002f) of a meridional return flow 
in the radiative zone in two other d ifferent models. The 
first o ne is just the same model of INandv fc Choudhuril 
l|2002j), but with the magneti c bouyancy mechanism of 
Dikp atTfc CharbonneavJ lll999l) . and is used mainly as 
reference considering both confined to convective zone 
and tachocline and deep meridional flow penetrating in 
the radiative zone. The second one introduces a bipoly- 
tropic approximation of the density profile to construct 
a meridional flow in the radiative and convective zones. 
This densi ty profile better resembles the sola r stan dard 
model of iBahcall. Pinsonneault fc Wasserbud lll995lk ?, 
and it is very interesting to explore how the hypotesis of a 
deeper meridional flow works there. In section 2 we show 
the mathematical formulation of the problem, based on 
the MHD induction equation, and define the meridional 
flow, differential rotation, magnetic diffusivity profiles and 
magnetig bouyancy mechanism to be used. Section 3 shows 
the results obtained with both models. Conclusions and 
a discussion of our results are exposed in section 4. Some 
details of the numerical implementation of the simulation 
are shown in the appendix A. 



where B(r,9,t) and V x (A(r,9,t)) correspond to the 
toroidal and the poloidal components of the magnetic field 
respectively; Q is the angular velocity, u = u r + ug is the 
velocity in the meridional plane and rj is the magnetic dif- 
fusivity. 

By replacing equations l|5J and in the induction 
equation Q and by separating the poloidal and toroidal 
components of the magnetic field, we obtain 

^ + i(u.V)(sA)=r ? (V 2 -^)A + S 1 (r,e,t) , (4) 

^ + ^§- r (™rB) + -?- e (ueB)] = (B p -V)n (5) 

-Vri x V x B + r?(V 2 - \)B , 

where s=r sin 6 and B P =V x A. It can be observed that 
a source term Si(r,8, B,/,) has been added (by hand) to the 
right side of eq. (QJ . This term is very necessary in our model 
for two reasons: it represents the magnetic buoyancy mech- 
anism which transports the toroidal field from the base of 
the solar convection zone to the surface, and it allows the 
surface regeneration of the poloidal magnetic field. We will 
discuss later about the functional form and physical content 
of this term. 



2.1 Differential rotation 

As we disscused before, the helioseismology gives a good 
characterization of the radial shear and the latitudinal dis- 
tribution of the solar angular velocity. It has been found a 
thin layer of substancial radial shear called tachocline, lo- 
cated at the base of the solar convection zone, where the 
toroidal magnetic field is generated (Q effect). An analyti- 
cal expression can be infer red from these mesurements, an d 
we will use the one used bv lDikpati fc Charbonneaul il999t) . 
who were the first to include a solar-like differential rotation 
profile in a kinematic model, 

Q(r,8) = n c + l[l + erf{2 T —- l£)](fi s (0) - fi c ) . (6) 
Z di 

Here, Q, s (8)=Q,Eq + a2 cos 2 # + a± cos 4 8 is the latitudinal 
differential rotation in the surface and erf(x) is an error 
function that confines the radial shear to a tachocline of 
thickness di=O.O5-R0. In this expression, a rigid core rotates 
uniformly with angular velocity f2 c /27r=432.8. Other values 
are f2 Eq /27r=460.7, a 2 /27r=-62.9, a 4 /27r=-67.13 nHz and 
r c =0.7R Q . 



2 MATHEMATICAL FORMALISM 

The MHD induct ion equation go verning the evolution of the 
magnetic field is JCowlinslll957h 

^ = Vx(UxB)+i)V 2 B . (1) 

By assuming spherical symmetry, the magnetic and velocity 
fields can be writen as 

B = B(r,8,t) + V x(A(r,8,t)), (2) 
U = u(r,0) +rsin0fi(r,0) , (3) 



2.2 Meridional circulation 

An analytical expression for the velocity in the meridional 
plane is more difficult to find that the previous one, because 
there is not enough observational data. On one hand, the 
observations suggest a poleward flow with an average ve- 
locity of 10-20 m/s dCiles et al.lll997D in the solar surface, 
but it is little what we know about the return flow, ex- 
cept that it must exist to satisfy mass conservation. On the 
other hand, numerical simulations of turbulent convection 
zones show that a structure of plumes in the base of the 
solar convection zone is able to push down the flow (and 
the magnetic field) to the radiative zone, where there is 
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a net movement towar ds the equator, with a mean veloc- 
ity of around 3 ms" _1 ferumell Hurlburt fc Toomrelll993 
iMiesch et al.ll200nF> . In this sense, we consider a single con- 
vection cell for each meridi onal quadrant, based on th e 
next equation introduced bvlDikpati fc Choudhuril lll994l) : 
IChoudhu ri 1 Schiisler fc Dikpatil lll995l) and previoslv used 
bv lNandv fc Choudhuril (120021) : 

(7) 



p(r)u= Vx [ip(r,e)e^] 
where ip is the stream function given by 

vr(r - R b ) 



iprsm8 = (r — Rb)ipo sin[ 
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/3 2 r(9-7r/2K 



(8) 



and p is the density profile for the sun. Most of the models 
use a profile for an adiabatic gaseous sphere with a specific 
heat ratio coefficient 7=5/3, which corresponds to a con- 
stant polytropic index m=l .5. So we have 



r 



0.95) r 



(9) 



where we chose C=3.60x 10 3 gr cm 3 as the surface density 
value dPinzon fc Calvo-Mozo|l200ll) . The coefficient ipo was 
chosen in such a way that the maximal latitudinal velocity 
at middle latitudes is 20 m s _1 . Other values in the eq. © 
are: /3i=1.65 x 10 10 cm" 1 , /3 2 =2.2xl0 10 cm" 1 , e=2. 0000001, 
r =(R e - R min )/4.15 and T=3.47 x 10 10 cm. Here R min is 
the minimal r coordinate value for the integration range, 
and the free parameter Rb is the max imal depth of return 
flow (See lDikpati fc Choudhuril lll994l) for more details). 



2.3 Magnetic buoyancy (MB) 

We mentioned above that a buoyancy magnetic force up- 
wards is generated in the base (or below) of the solar con- 
vection zone, due to of the density difference between the 
inner and the outer regions of a magnetic flux tube. Simula- 
tions in this sense have shown that toroidal magnetic tubes 
of 10 G emerge to the surface and are twisted by Corilis 
forces to form tilted bipolar active magnetic regions (tilted 
BMRs). They have also shown that the tilt and the latitud 
of those BMRs strongly depend of the initial value of the 
toroidal magnetic field. Initial values greater than 1.6 x 10 5 
G don't produce tilt in BMR, and values lower than 6 x 10 4 
emerge radially with t ilts in disagreement with Joy's law 
JP'Silva fc Choudhuril Il993t iFan. Fisher fc DeLucal Il993t 
ICaligari. Moreno-Insertis fc Schiisleilll995lll99^ . We intro- 
duce a simplified form of inc luding these results, due to 
iDikpati fc Charbonneaul il999Tl . with an initial value of 10 5 
G. Its analytical form is 



Si(r,B;B) 



^BMm + erfi^p-] 



(10) 



x [1 - er/( 



r - r 3 , 
d 3 ■ 



x [i + ( £%M )2] - lcos0 . 

Bo 

It can be seen that this term acts non-locally in B 
<Durnevll995l. Il996l . ll99H) : values of toroidal magnetic field 
B at the tachocline produce proportional poloidal magnetic 



fields at a thin layer close to the surface, defined by the error 

anihillates 



3{r c> B,t) -,2!-! 



(erf) functions. The last term, [l + (- 
the subsequent increase of A at the surface, beyond some 
maximal level Bo- It is, in other words, a saturation term, 
and it is the only source of non-linearity in the system. The 
other parameters are r 2 =O.957?0, r 3 =R Q , d2=d3=0.025-R Q , 
so the alpha effect is confined beneath the surface in a layer 
of thickness 0.025.R©. 



2.4 Magnetic diffusivity 

Magnetic diffusivity r\ is different for the radiative and 
the convective zones. The turbulent regime present at the 
convective zone makes r\ two ord ers of magnitude higher 
than t he radiative one (see fig 1 of lDikpati fc Charbonneaul 
J1999D 1. 

v(r)=Vc+ 1 f[l + erf(2 r -^)] , (11) 
with r/ c =2.2 x 10 9 and ^=0.5x10" cm 2 s" 1 and r c =0.7R Q . 



3 RESULTS 

We solved the model described above by means of the ADI 
method on a two-dimensional mesh of 64x64 spatial divi- 
tions, with 0.557?© ^ r ^ IRq and 0^8^ 7r/2, and by 
imposing the following boundary conditions: 



at 
at 



= 



W9- 9A 



at r = 0.55.Rq 
at r = Rq : 



A = 0, B = , 
0, B = 

A = 0, B = , 
B = . 



(12) 
(13) 
(14) 
(15) 



At the upper boundary it must be ensured, in addition, that 
A satisfies the free space condition, 



(V 2 



-)A = 



(16) 



The initial condition is not determinant for the evolution of 
the model, as long as the So coefficient is large enough to 
avoid the magnetic field to diffuse. In this case the system 
always relaxes to peridic solutions. In all cases, we took the 
initial condition 



B(r,6,0) = Bosin(26>) 



(17) 



T he code was tested by using al l profiles and parameters 
of lDikpati fc Charbonneaul l|l999F) . and by reproducing the 
results therein. The numerical details on discretization and 
integration are shown in appendix A. 

3.1 Meridional flow confined to the convective 
zone 

In this first model, the meridional flow penetrates until 
0.675-Rq, that is, it is confined to the solar convection zone 
and the tachocline. The free parameters used in this model 
are shown on the left side of Tabled As the butterfly dia- 
gram shows (Fig^V), the maximal toroidal magnetic field 
in the base of SCZ is located at latitudes between 75 and 
90 degrees, with weak branches migrating in the equator- 
ward direction. Under the primordial assumption that the 
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Table 1. Free parameter values for a meridional flow confined to 
the convective zone 



Parameter Value Result Value 

U 2000 cm s" 1 B^max 2.83 x 10 5 G 

So 5 cm s — 1 B r niax 150 G 

r\ T 0.8 X 10 11 cm 2 s" 1 T 27.8 years 
R b 0.675R© 



Table 2 . Free parameter values and results for a deep meridional 
flow 



Parameter Value Result Value 

U 2000 cm s" 1 B^max 2.58 x 10 5 G 

So 5 cm s — 1 B r max 50 G 

r\ T 0.8 X 10 11 cm 2 s" 1 T 28.8 years 

R b 0.61R Q 
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Figure 1. Butterfly diagrams for the, A, toroidal magnetic field 
at the base of solar convection zone (r=0.7i?Q) and, B, radial 
field at the surface (r=iio), by using the parameters shown in 
Table The contours are equally spacied, with solid (dashed) 
lines for positive (negative) values. Time in years and latitude in 
degrees. 



toroidal field lines emerge to the surface like magnetic flux 
tubes, this model generates sunspots close to the poles. The 
rigth side of Table Q summarizes the most important results 
of this model. We obtain a period slightly larger than the 
one observed and the maximal intensities of the magnetic 
fields are in the correct order of magnitude for the toroidal 
one but one order the magnitude over the one expected for 
the radial field. 



3.2 Deep meridional flow 

A different situation is observed if the meridional flow pene- 
trates into the radiative zone, as deep as r=O.61i?0. A better 
distribution of the toroidal lines is obtained, however, a po- 
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Figure 2. Butterfly diagrams for the, A, toroidal magnetic field 
in the base of solar convection zone (r=0.7-Ro) and, B, radial 
field in the surface (t=Rq) with the parameters showed in table 
[5J Time in years and latitude in degrees. 



lar br anch still remains, in contrast to lNandv fc Choudhuril 
( 2002). There are two zones of maximal intensity: one below 
45 degrees and another above 80 degrees (at the poles) . This 
result can be interpreted in the following way: if most of the 
magnetic field is dragged below the tachocline, only a small 
portion emerges to the surface by the action of the magnetic 
bouyancy force. The rest migrates to the equator, dragged 
by the solar plasma, and will undergo the alpha effect only 
at middle and low latitudes. 

It is important to notice that this model gives a period 
of 28.8 years (close to the observed one) and maximal val- 
ues of both toroidal and radial fields within the right order 
of magnitudes. The parameters we used and the results are 
summarized in Table [5] The large scale behaviour of the 
model can be observed in the butterfly diagrams of Fig 2. 
The equatorward migration of the branches of toroidal mag- 
netic field is evident. 
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Table 3. Free parameter values and results for a more realistic 
density profile model 



Parameter 


Value 


Result 


Value 


U 


2000 cm s" 1 




3.77 x 10 5 G 


So 


5 cm s — 1 


_B r max 


300 G 


r)r 


0.8 X 10 11 cm 2 s" 1 


T 


72.1 years 




0.61.Rq 








0.6 0.7 0.8 0.9 1 
r 

Figure 3. Density profiles in the solar interior for the usual way 
(dashed line) and the bipolytropic model (solid line). In both 
graphics the density has been normalized to a surface value of 
3.603 X 10~ 3 g cm -3 like the one in the bipolytropic profile. Den- 
sity in g cm -3 , r in Rq. 

3.3 A more realistic density profile 

As the model of Nan dy and Choudhouri does 
iNandv fc Choudhuril 120021) . the model above employs 
a single polytropic density profile (Eq. [5J|, with m=1.5 
(7=5/3). This is just true for the convective zone (r>0.7), 
but not for the radiative one. We can introduce a more 
realistic densit y profile by using a bipoly tropic aproxi- 
mation due to IPinzon fc Calvo-Mozol i200ll) . with 7=5/3 
for the convective zone and 7=1.26 for both radiative 
zone and core (Fig El. We fit the numerical results of 
IPinzon fc Calvo-Mozol (1200 1*1 by the analytical expression 

PW = (— - 0.904) 2 ' 4 , (18) 
r 

with the same surface density value as 

Despite that this profile closely resembles the 
internal structure of the solar stan dard mode l 

of iBahcall. Pinsonneault fc Wasserbuel lll995l) : 

iBahcall.Phisorineault^^fj^sul I2OO1T I our results are 
not so good as before. 

On one hand, we obtain again two maximal intensity 
zones: the first one very close to the poles, but the second 
one at middle latitudes, between 60 and 40 degrees, in dis- 
agreement with the observations. On the other hand, since 
the bipolytropic density profile is around three times larger 
that the usual density profile at the return point, r=O.617?0, 
(see Fig[3J, the counter flow is around three times lower 
that before in order to hold mass conservation, V • (pU)=0. 
Thus, our model gives a new period much larger than ob- 
served (actually, three times larger, but this is not a linear 
relationship). If we allow the flow to go deeper, the second 
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Figure 4. Butterfly diagrams for the toroidal magnetic field in 
the base of solar convection zone (r=O.7i?0) for a meridional flow 
penetrating into the radiative zone and a bipolytropic density 
profile. 

zone of maximum toroidal field approaches the equator, but 
at the cost of increasing the period very quickly. 



4 CONCLUSIONS 

In this work three different models of solar dy- 
namo are shown. We dev eloped these models by us - 
ing the velocity field of iNandv fc Choudhuril <2002l) . 

the magnetic buoyan cy mechanism introduced by 

Dikp ati fc Charbonneaul l)l999h and different values of 
the diffusivity coefficie nt for the convective a nd th e 
radiative zones, as in [ Dikpati fc Charbonneaul £l999); 
INandv fc Choudhuril feOOaT In other words we have recov- 
ered the TNEmd^^^Choudhuril d2002l) results even with a 
different source formulation. This suggests that the good 
performace of the kinematic models in the papers above is 
general in nature. Moreover, our results suggests that, at 
least for this kind of models, the effect of a deep meridional 
flow to solve the apperance of sunspots at high latitudes is 
general, too. 

The most important results are: 

(i) If the meridional flow is confined to the convective 
zone (r>O.675i?0), the emergence latitude of sunspots is 
just near the poles, as expected. This model also gives wrong 
maximal values of the radial field intensities. 

(ii) When the flow penetrates until O.6I-R0, the model 
gives a better solution. However, two regions of maximal 
toroidal fi eld appear (in contrast with the single region ob- 
tained bv INandv fc Choudhuril J2002|) 'l: one of them is lo- 
cated in the right range of latitudes and the other remains 
very near to the poles. Except for this last region, the model 
approximates well the observed butterfly diagram. 

(iii) In the third model we attempt to improve it by in- 
corporating a more realistic density profile, but the results 
are not satisfactory. There are two zones of maximal toroidal 
field as before, but the second one appears at latitudes to 
high to agree with observations. Moreover, the period of cy- 
cle is three times larger than observed. 

Our simulations show again the good effects of a deep 
meridional flow to solve the appearance of sunspots at 
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high latitudes, proposed by iNandv fe Choudhuril (§002), 
but this hypothesis deserves a more extended discussion. 
First, recent results on magn etic tube flux simulations 
iRempel. Schiiser fe Totbl l200dl suggests that a subadia- 
batic stratification may suppress the magnetic buoyancy 
force and address the equilibrium of a to roidal band. De- 
spite that iRempel. Schiiser fc Totbl i2000f) do not speak of 
a deep flow, such kind of stratification is the one to be found 
in the radiative zone and the tachocline, supporting the 



process suggested by INandv fe Choudhuril (l2002f) to drag 
the toridal field towards lower latitudes. However, they do 
not use a subadiabatic stratification, but an adiabatic back- 
ground (7=5/3) everywhere. We have introduced a subadi- 
abatic background in a simple way, by using a more realistic 
solar structure with a bipolytropic density profile, but our 
results are not in agreement with the observations. Never- 
theless, this may suggest that a more realistic meridional 
circulation profile is needed, not that the idea of a mag- 
netic buoyancy suppressed by a sub-adiabatic stratification 
should be rejected. 

Second, it should be noticed that the discussion around 
a deep meridional flow hipothesis is not closed. On one 
hand, turbulent convection models obtain flow penetra- 
tion (in a plumes dominated structure ) in the stable layer 
l|Miesch et al.ll2000t Irbbias et alJl200lT l. and some analysis 
on helioseismological observations give a pole ward merid- 
ional flow across the entire convective zone llGiles et alJ 
Il997f) . supporting the idea of a deep flow. On the other hand, 
the meridional flow is such a weak flow that it is very unlikely 
that it can penetrate the strongly subadiabatically strati- 
fied radiative core of the Sun. In addition, other problems 
would have to be reconsidered if a deep flow is assumed, like 
a larger angul ar momentum transfer to the radiative zone 
iDurnevl2000l) and the changes in the relative abundance of 
Lythium a nd other elements in both radiative and convec - 
tive zones dZahr]|200l|: iBrun. Turck-Chieze fc Zahnlll999l) . 
but these are still open problems. Despite these uncertain- 
ties, our aim with this paper is to verify what would happen 
if this scenario were possible. 
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APPENDIX A: NUMERICAL METHODS 

The equations and ^ are two coupled partial differen- 
tial equations of advection-diffusion type, with a non-linear 
term (Si). There are two important aspects to have into ac- 
count: numerical dicretization and numerical integration. In 
the first one, the advective-diffusive character of the equa- 
tions must be considered. In order to ensure precision and 
stability for all times the advective terms are discretized 
by using the Lax-Wendroff mechanisms, which is a second- 
order discretization in both space and time, and the diffusive 
terms are discretized by a simple forward-time centered - 
space (FTCS) mechanism jAmeslll977t iPress et allll992Fl . 
In the second one, the two dimmensional character and the 
non-linearity are considered. Alternating Direction-Implicit 
(ADI) met hod is a good numerical techni que for these char- 
acteristics dAmeslll977i IPress et al]ll992l) . 

A first step is to write the equations in operator nota- 
tion, as follows: 



dA 

— = (L r + Le)A + Si , 
^ = (L r + L e )B + S 2 , 
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and Si and S2 are the crossed terms. Si is given by 
equation 1 1 lit and S2 is given by 

r *sm0(B p .V)n = -— sin0(-A-r — ) + — (Acot0+ — ). 

(A7) 

In the ADI method, the time step is divided into two 
steps of size At/2. In each half step, one spatial dimmen- 
sion is treated implicitly and the other is treated explicitly. 
We treated the 9 terms in implicit form and the r terms in 
explicit form in the first half step, as follows: 



ij 2 = L e A n+1/2 + L r A n + S1", (A8) 



£,n + l/2 _ j^n 

13 M/2 lj = L e B n+1/2 + L r B n + S> 2 \ (A9) 

where i,j = 0, 1, N are the spatial divisions and n is 
the time step. The next half step is treated in the reverse 
manner, 

A n+1 — A n+X l 2 
13 At/ ^ J = L e A n+1/2 +L r A n+1 + S" +1/2 , (A10) 

nn+l _ n"+l/2 

13 <j = L e B n+1/2 +L r B n+1 + S 2 n+1/2 . .(All) 

The source terms are always treated explicitly, to guar- 
antee the linearity of equations in A and B. These equations 
can be organized in such a way that they can be solved 
by a standard tridiagonal algorithm. The numerical treat- 
ment of the boundary cond itions is developed in detail in 
iDikpati fc Choudhuril lll995|) . 
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